MODIFIED SHALLOW WATER EQUATIONS FOR SIGNIFICANTLY 

VARYING BOTTOMS 



DENYS DUTYKH* AND DIDIER CLAMOND 

Abstract. In the present study we propose an modified version of the nonlinear shahow 
water (Saint- Venant) equations for the case when the bottom undergoes some significant 
variations in space and time. The model is derived from a variational principle by choosing 
the appropriate shallow water ansatz and imposing some constraints. Our derivation pro- 
cedure does not explicitly involve any small parameter and is straightforward. The novel 
system is a non-dispersive, and non-hydrostatic extension of the classical Saint- Venant 
equations. We also propose a finite volume discretization of the obtained hyperbolic sys- 
tem. Several test-cases are presented to highlight the added value of the new model. Some 
■ implications to tsunami wave modelling are also discussed. 
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1. Introduction 

The celebrated classical nonlinear shallow water or Saint- Venant equations were derived 
for the first time in 1871 by A.J.C. de Saint-Venant [17], an engineer working at Ecole 
Nationale des Fonts et Chaussees in France. Currently these equations are widely used in 
practice and the literature counts many thousands of publications devoted to the applica- 
tions, validations and numerical solutions of these equations [2, 34, 35, 38, 75, 88]. 

Some important attempts have been also made to improve this model from physical 
point of view. The main attention was payed to various dispersive extensions of shallow 
water equations. The inclusion of dispersive effects resulted in a big family of the so-called 
Boussinesq-type equations [21, 26, 32, 54, 57, 58, 60, 63]. Many other families of dispersive 
wave equations have been proposed as well [14, 18, 20, 42, 56, 67, 81]. 

However, the interaction of the waves with mild or tough bottoms has always attracted 
the particular attention of researchers [58, 6, 13, 33]. There are a few studies which attempt 
to include the bottom curvature effect into the classical Saint-Venant [17, 72] or Savage- 
Hutter^ [66, 41, 86] equations. One of the first studies in this direction is perhaps due 
to Dressler [22]. Much later, this research was pursued almost in the same time by 
Keller, Bouchut and their collaborators [49, 11]. We note that all these authors used 
some variants of the asymptotic expansion method. The present study is a further attempt 
to improve the classical Saint-Venant equations by including a better representation of 
the bottom shape. Moreover, as a derivation procedure we choose a variational approach 
based on the relaxed Lagrangian principle [16]. The derivation presented below was already 
communicated by the same authors in a short note announcing the main results [25]. In 
this study we investigate deeper the properties of the proposed system along with its 
solutions through analytical and numerical methods. The obtained results may improve 
the modeling of various types of shallow water flows such as open channel hydraulics [15], 
atmospheric and oceanic flows [61], waves in coastal areas [9, 5, 30, 32], and even dense 
snow avalanches [66, 47, 86]. 

The present article is organized as follows. After some introductory remarks, the paper 
begins with the derivation and discussion of some properties of the modified Saint-Venant 
(mSV) equations in Section 2. Then, we investigate the hyperbolic structure and present 
a finite volume scheme in Section 3. Several numerical results are shown in Section 4. 
Finally, some main conclusions are outlined in the last Section 5. 

2. Mathematical model 

Consider an ideal incompressible fluid of constant density p. The horizontal independent 
variables are denoted by a; = {xi^X2) and the upward vertical one by y. The origin of the 
Cartesian coordinate system is chosen such that the surface |/ = corresponds to the still 
water level. The fluid is bounded below by the bottom aXy = —d{x, t) and above by the free 
surface aX y = ri{x,t). Usually, we assume that the total depth h{x,t) = d{x,t) + ri{x,t) 



The Savage-Hutter equations are usually posed on inclined planes and they are used to model various 
gravity driven currents, such as snow avalanches [4]. 
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Figure 1. Definition sketch. 



remains positive h{x, t) ^ /iq > at all times t G [0, T]. The sketch of the physical domain 
Q X [— d, f2 C R2 is shown on Figure 1. 

Traditionally in water wave modelling the assumption of flow irrotationality is also 
adopted. Under these constitutive hypotheses, the governing equations of the classical 
water wave problem are [52, 73, 55, 84, 48]: 



+ dy'<P 

dtv + (V0) • (Vr^) - dy(P 

dtd + {yd) ■ + dy(f) 



0, {x,y)enx[-d,r]], (2.1) 

0, y = v{x,t), (2.2) 

0, y = v{x,t), (2.3) 

0, y = -d{x,t), (2.4) 



where is the velocity potential (by definition, u = V0 and v = dy(j)), g > is the 
acceleration due to gravity and V = {dx^, 0^2) denotes the gradient operator in horizontal 
Cartesian coordinates. 

The assumptions of fluid incompressibility and flow irrotationality lead to the Laplace 
equation (2.1) for the velocity potential (j){x,y,t). The main difficulty of the water wave 
problem lies on the boundary conditions. Equations (2.2) and (2.4) express the free-surface 
kinematic condition and bottom impermeability, respectively, while the dynamic condition 
(2.3) expresses the free surface isobarity. 

It is well-known that the water wave problem (2.1) - (2.4) possesses several variational 
structures [64, 83, 53, 87, 12]. Recently, we proposed a relaxed Lagrangian variational 
principle which allows much more freedom in constructing approximations compared to 
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the classical formulations. Namely, the water wave equations can be obtained as Euler- 
Lagrange equations of the functional jjj ^d^xdt involving the Lagrangian density [16]: 

^ = {dtri + ii-Vrj-v)^ + {dtd + fi ■ V d + v) '(j) - \g7f 

/V 
[^^^u - ^u^ + uv - + (V • /X + 9yZ/)0] d?/, (2.5) 

where the over 'tildes' and 'wedges' denote, respectively, a quantity traces computed at 
the free surface y = ri{x,t) and at the bottom y = —d{x,t) (we shall also denote below 
with 'bars' the quantities averaged over the water depth); {u, v, ^t, u} being the horizontal 
velocity, vertical velocity and associated Lagrange multipliers, respectively. The last two 
additional variables {/x, z/} are called the pseudo- velocities. They formally arise as Lagrange 
multipliers associated to the constraints u = V0, v = (py However, once these variables 
are introduced, the ansatz can be chosen regardless their initial definition. The advantage 
of the relaxed variational principle (2.5) consists in the extra freedom for constructing 
approximations. 

2.1. Constrained shallow water ansatz. In order to simplify the full water wave prob- 
lem, we choose some approximate but physically relevant representations of all dependent 
variables. In this study, we choose a simple shallow water ansatz, which is a velocity field 
and velocity potential independent of the vertical coordinate y such that 

^ (f){x,t), u = ^ u{x,t), V = u ^ v{x,t), (2.6) 

where v{x,t) is the vertical velocity at the bottom. In this ansatz, we take for simplicity 
the pseudo-velocities to be equal to the velocity field u = fi, v = u. However, in 
other situations they can differ (see [16] for more examples). The Lagrangian density (2.5) 
becomes: 

^ = {dth + u-Vh + hV -u)^ - \grf + \h{u^ + v"^), (2.7) 
where we make appear the total water depth h = rj + d. 

Remark 1. Note that for ansatz (2.6) the horizontal vorticity uj and the vertical one ( 
are given by: 

Consequently the flow is not exactly irrotational in general. It will be confirmed below one 
more time when we establish the connection between u and V0. 

Now we are going to impose one constraint by choosing a particular representation of 
the fluid vertical velocity v{x, t) at the bottom. Namely, we require fluid particles to follow 
the bottom profile: 

V = -dtd - u-Vd. (2.8) 

The last identity can be also seen as the bottom impermeabiliy condition within ansatz 
(2.6). After substituting relation (2.8) into Lagrangian density (2.7), the minimization 
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procedure yields: 

50 : = dth + V ■ [hu], (2.9) 

Su : = u - V(f) - vVd, (2.10) 

Sri : = dt(l) + gri + u-V(f)- l{u^ + v^). (2.11) 

Taking the gradient of (2.11) and eliminating of (f) from (2.10) gives us this system of 
governing equations: 

dth + V ■ [hu] = 0, (2.12) 
dt[u-vVd] + V[gri + + \v'^ + vdtd\ = 0, (2.13) 

together with relations: 

u = V0 + vVa, V = — Otd — u-Vd = — . 

1 + |Vap 

Remark 2. The classical nonlinear shallow water or Saint- Venant equations [17, 72] can 
be recovered by substituting v = into the last system: 

dth + V -[hu] = 0, 

dtu + V[gri + | w^] = 0. 

The last equation is generally written in the literature into the conservative momentum flux 
form [19, 65, 34].- 

dt[hu] + V[hu^ + Igh^] = ghVd. 

2.2. Properties of the new model. From governing equations (2.12), (2.13) one can 
derive an equation for the horizontal velocity u and for the momentum flux hu: 

dtU + u-Vu + gVr] = 'jVd + u AVv AVd, (2.14) 

dt[hu] + V[hu^ + \gh^] = {g + -f)hVd + huAVvAVd, (2.15) 
where 7 is the vertical acceleration at the bottom defined as: 

= ^ = dtv + {u-V)i!. (2.16) 

Remark 3. Note that in equations (2.14) and (2.15) the last term on the right-hand 
side cancels out in one horizontal dimension. It can be seen from the following analytical 
representation which degenerates to zero in one horizontal dimension: 

hu AVv AVd = (Vv) {hu -Vd) - (Vd) {hu -Vv). 

This property has an interesting geometrical interpretation since hu A Vt; A Vc? is a hori- 
zontal vector orthogonal to u and thus vanishing for two-dimensional waves. 

One can also derive an equation for the energy flux: 



InP + -u^ r7^ — d"^ 
^ +9 ^ 



+ V 



tip + v"^ 



+ gr] I hu 



+ -f)hdtd. (2.17) 
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Obviously, the source term on the right-hand side vanishes if the bottom is fixed d = d{x) 
or, equivalently, dtd = 0. 

The last conservation law is closely related to the Hamiltonian structure of the mSV 
equations (2.9)-(2.11). Namely, these equations possess a Hamiltonian structure with 
canonical variables h and 0, i.e., 

d]i _ 6H 50 _ _6n 
'dt ~ ~5^' 'dt ~ ~lh' 

where the Hamiltonian is 

1 /" r ,n9 ,,-,t,9 hldtd + V(f)-Vd]'^'] 

^ = ly \^9{h-d)' + h\V<j>\' - ^^^—^^^^^1 d'x. (2.18) 

One can easily check, after computing the variations, that the Hamiltonian (2.18) yields 

dfh = -V 



dt^ = -g{h~d) - i|V^ 



,2 [dtd + V(t)-Vd]'^ 



l + |Vd|2 

which are equivalent to the system (2.9) - (2.11). 

Remark 4. If we rewrite Hamiltonian (2.18) in the following equivalent form: 

n = I j [gr]^ + hu^ + h{v + dtdf - h{dtdf] d?x, (2.19) 

one can see that Hamiltonian (2.18) is actually positive definite if the bottom is static, i.e., 
d = d{x) or dtd = 0. In other words if there is no external input of energy into the system. 

Moreover, the Hamiltonian structure of the classical Saint- Venant equations can be re- 
covered substituting v = into the last Hamiltonian (2.19).- 

Ho = ^ j {gV^ + hu^} d^x. 

2.3. Steady solutions. We consider here the two-dimensional case in the perspective to 
derive a closed form solution for a steady state fiow over a general bathymetry. We assume 
the following upstream conditions at x — )■ — oo: 

?7 — 7- 0, d ^ do, u — 7- uo ^ 0. 

Physically, these conditions mean that far upstream we consider a uniform current over a 
horizontal bottom. The mass conservation in steady condition yields 

hu = doUo, 

while the momentum conservation equation becomes 

gh + lu^[l+ {d,df] = gdo + ^u^ 

Combining the last two relations yields the following cubic equations for the total water 
depth 

gh' - (gdo + U^)h^ + i(rfoWo)'[(l + {d^df] = 0. (2.20) 
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As it is well known, a cubic equation has three (in general complex) roots. From physical 
interpretation of h, we are only interested in real positive roots. 

Remark 5. It is straightforward to derive a similar equation for steady solutions to the 
classical Saint- Venant equations 

gh^ - {gdo + |mo)/i^ + |(rfoMo)^ = 0. 

The last relation can be also obtained from equation (2.20) by taking d^d = 0. Consequently, 
we can say that steady solutions to classical Saint- Venant equations do not take into account 
the bottom slope local variations. 

One solution of (2.20) is 

h = do/[AC-^/^ + |C^/^] , C = 3 - 3^3 - 3B, (2.21) 



with 



1 + (5.c/)2 ^ ' ' ' 1 + {d^dy 

This solution is real only if ^ 3^^ (in that case equation (2.20) has only one real root) 
or equivalently if 



27F^djf - (f2-1)2(f2 + 8) ^0, F ^ % / /MT, 

where F is a Froude number. In the case of the classical Saint- Venant equations (i.e., taking 
dxd = 0), the last condition is satisfied only for the critical regime F = 1. In the modified 
model solution, (2.21) is real in a vicinity of F = 1, which depends on the magnitude of 
the bottom local slope dxd. However, if this solution is real, it is necessary negative — 
see the definition of C in (2.21). Hence, it is unphysical and must be rejected. In other 
words it means that under condition B^ > 3A^, the flow cannot be steady. Physically, a 
singularity such as wave breaking does not allow the flow to reach the steady state. 

If B^ < 3A^, the equation (2.20) has three real roots, two of them being positive. The 
first root h = corresponds to the subcritical regime: 



/i"*" = do 



2 v^cos(| arccos (-3-i/2^yl-3/2) - |7r) 



and one supercritical root h = h 



h = dn 



2 cos (I arccos {-3-^^^BA-^/^)) 



We note that h~ < h~^. 

Finally, in the limiting case B^ = 3A^, for a given Froude number F , there is only one 
absolute value of the slope for which this identity is satisfied, that is 

(dxdy = (F2-1)2(f2 + 8)/(26F2). (2.22) 

For instance, if d^d = then B^ = 3A^ if and only if F =1. 
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Gravity acceleration g: 


1 m s ^ 


Undisturbed water depth (Iq: 


Im 


Deformation amplitude a: 


0.5 m 


Half-length of the uplift area h: 


2.5 m 


Upstream flow speed, Uq 


2.0 ms-^ 



Table 1. Values of various parameters used for the steady state computation. 

Consider now the opposite case when the slope is given 7^ 0. In this case condition 
(2.22) is satisfied for two values of the Froude number F = F^: 

{F-f = 6 + (9,rf)2cos(|arccos (-[1 + [d^df]-^/^) - fvr) -2^1, 
= 6 + (9,rf)2cos(|arccos (-[1 + [d^df]-^'^)) -2^1. 

IfF~<F <F + , we have > 3A'^ and there are no steady states. Therefore, we must 
have F ^ F ^ or F 

In order to illustrate the developments made above, we compute a steady flow over a 
bump. The bottom takes the form 

d{x) = do - ab-^{x^ - b^fR{b^-x^), 

where H(x) is the Heaviside step function [1], a and b being the bump amplitude and its half- 
length, respectively. The values of various parameters are given in Table 1. We consider 
here for illustrative purposes the supercritical case for the classical and new models. The 
result are shown on Figure 2 where some small differences can be noted with respect to 
the classical Saint- Venant equations. 

3. Numerical methods 

In this Section we discuss some properties of the system (2.12), (2.13) and then, we 
propose a space discretization procedure based on the finite volume method along with a 
high-order adaptive time stepping. 

3.1. Hyperbolic structure. From now on, we consider equations (2.12), (2.13) posed in 
one horizontal space dimension for simplicity: 

dth + d^[hu] =0, (3.1) 
dt[u - vd^d] + [gr] + + \v'' + vdtd\ = 0. (3.2) 

In order to make appear conservative variables we will introduce the potential velocity 
variable U = d^ct)- From equation (2.10) it is straightforward to see that U satisfies the 
relation 

U = u — V d^d, 
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Figure 2. Supercritical steady state solutions over a bump for the Froude 
number F =2. Comparison between the classical and modified Saint- Venant 
equations. 



Depth averaged and vertical bottom velocities can be also easily expressed in terms of the 
potential velocity U 



u 



U - jdtd) {d^d) 

1 + {d^dy ■ 



dtd + U d^d 

' 1 + {d^dy ■ 



Consequently, using this new variable equations (3.1), (3.2) can be rewritten as a system 
of conservation laws 



dth + 



h 



U - (dtd) {d,d) 



dtU + 



g{h-d) + 



1 + (d^dy 

lU^ - 2U{dtd) {d^d)-{dtdf 

2 1 + {djy 



0, 
0. 



For the sake of simplicity, we rewrite the above system in the following quasilinear vectorial 
form: 



dtw + d^f{w) = 0, 



(3.3) 
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where we introduced the vector of conservative variables w and the advective flux f{w): 



( 



w 



/H 



h 



g{h -d) + 



U-{dtd){d^J) \ 

1 + {d.,df 
U^-2U{d,d){dtd)-{dtdY 



, 2 [1 + id^d^] 

The Jacobian matrix of the advective flux f{w) can be easily computed: 



A{w) 



df{w) _ 1 \U~{dtd){dJ) h 
dw ~ l + {d,d)^[9{l + {d,dy) U-{dtd){dJ) 



u 



1 + {dJY 
g u 



The matrix k.[w) has two distinct eigenvalues: 

U - [dtd) {d^d) 



A 



± 



± c 



M ± C, 



gh 



Remark 6. Physically, the quantity c represents the phase celerity of long gravity waves. 
In the framework of the Saint- Venant equations, it is well known that c = y/gh. Both 



expressions differ by the factor a/1 + {d^dy. In our model, the long waves are slowed 
down by strong bathymetric variations since fluid particles are constrained to follow the 
seabed. 

Right and left eigenvectors coincide with those of the Saint- Venant equations and they 
are given by the following matrices 



R 



-h 
y/gh 



h 

'gh 



(gh)-'" 



Columns of the matrix R constitute eigenvectors corresponding to eigenvalues A~ and 
respectively. Corresponding left eigenvectors are conventionally written in lines of the 
matrix L. 

3.2. Group velocity. We would like to compute also the group velocity in the framework 
of the modified Saint- Venant equations. This quantity is traditionally associated to the 
wave energy propagation speed [73, 28]. Recall, that in the classical linearized shallow 
water theory, the phase c and group Cg velocities are equal [73]: 

OJ r-^ dw 



gh, 



-JT = ygh 



k ' ' dk 

where u = ky/gh is the dispersion relation for linear long waves, k being the wavenumber 
and uj being the angular frequency. 

In order to assess the wave energy propagation speed we will consider a quasilinear 
system of equations composed of mass and energy conservation laws: 

U - {d,d){dtd)' 



dth + 



h 



dtE + 



{E + \gh^ 



1 + {dJY 
U-{d,d){dtd) 
1 + {d.df 



0, 



-{g + l)hdtd, 
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where 7 is defined in (2.16) and E is the total energy considered aheady above (2.17): 
u^ + v^ g jr]^ - d^) _ h U^ + {dtdf g {h^ - 2hd) 
2 ^ 2 ~ 2 1 + {d^df ^ 2 

The last formula can be inverted to express the potential velocity in terms of the wave 
energy: 

= [1 + [d^df] -gh + 2gd^ - {dtdf. 

In the spirit of computations performed in the previous section, we compute the Jacobian 
matrix J of the mass-energy advection operator: 



1 



U - {d^d){dtd) + h^^ 



dh BE 
gh[U - (d^dMd)] + {E + hh^)^ U - {d,d){dtd) + {E + \gh^^ 



1 + {d^dy 

where partial derivatives are given here: 



dU , gh^ + 2E dU l + id^df 

[1 + [d^d) ] 



dh 2h'^U ' dE hU 

Computation of the Jacobian J eigenvalues leads the following expression for the group 
velocity of the modified Saint- Venant equations: 

2 _ gh U - {d^d){dtd) 

~ 1 + {d^dy U 

The last formula is very interesting. It means that in the moving bottom case, the group 
velocity Cg is modified and does not coincide anymore with the phase velocity = gh[l + 
{dxdY]~^ . This fact represents another new and non-classical feature of the modified Saint- 
Venant equations. The relative difference between phase and group velocities squared is 

C2 U ' ^ ' 

which is not necessarily always positive. When it is negative, the energy is injected into 
the system at a higher rate than can be spreaded, thus leading to enenergy accumulation 
and possibly favouring the breaking events. 

3.3. Finite volume scheme. Let us fix a partition of R into cells (or finite volumes) 
Ci = , with cell center Xi = |(a;j_i -|- x^^^^i), i G Z. Let Axi denotes the length 
of the cell C,. Without any loss of generality we assume the partition to be uniform, i.e., 
Axi = Ax, Wi G Z. The solution w{x,t) is approximated by discrete values and, in order 
to do so, we introduce the cell average of w over the cell Ci, i.e., 

"^iW = TT / w{x,t)dx. 
Ah Jc^ 

A simple integration of (3.3) over the cell Ci leads the following exact relation 

dwi f {w{x,^i_,t)^ - f {w{Xi_,_,t)^ 
"dT ^ Ax 
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Since the discrete solution is discontinuous at cell the interfaces z G Z, the heart 

of the matter in the finite volume method is to replace the flux through cell faces by the 
so-called numerical flux function 



where w^fl are reconstructions of conservative variables w from left and right sides of 

each cell interface [8, 80]. The reconstruction procedure employed in the present study is 
described below. 

In order to discretise the advective flux f{w), we use the so-called FVCF scheme [39, 40] 

f{v) + fjw) , /H - fjv) 
J^(v,w) = S{v,w) . 

The first part of the numerical flux J^{v, w) is centred, while the second part is the up- 
winding introduced according to local waves propagation directions 

S{v,w) = sign ' sign(A) = i? • diag(s~, s"*") • L, = sign(A^). 

After some simple algebraic computations one can find expressions for sign matrix S coef- 
ficients 



is' -s ')^~gjh \s' + s 



all coefficients being evaluated at the average state of left and right face values. 

Taking into account the developments presented above, the semi-discrete scheme takes 
the form 



dt Ax ^ ' 

The discretisation in time of the last system of ODEs is discussed in Section 3.5. Meanwhile, 

we present the employed reconstruction procedure. 

3.4. High-order reconstruction. In order to obtain a higher-order scheme in space, 
we need to replace the piecewise constant data by a piecewise polynomial representation. 
This goal is achieved by various so-called reconstruction procedures, such as MUSCL TVD 
[51, 79, 80], UNO [46], ENO [45], WENO [85] and many others. In our previous study on 
Boussinesq-type equations [32], the UN02 scheme showed a good performance with low 
dissipation in realistic propagation and runup simulations. Consequently, we retain this 
scheme for the discretisation of the modified Saint- Venant equations. 

Remark 7. In TVD schemes, the numerical operator is required (by definition) not to 
increase the total variation of the numerical solution at each time step. It follows that the 
value of an isolated maximum may only decrease in time which is not a good property for 
the simulation of coherent structures such as solitary waves. The non-oscillatory UN02 
scheme, employed in our study, is only required to diminish the number of local extrema in 
the numerical solution. Unlike TVD schemes, UNO schemes are not constrained to damp 
the values of each local extremum at every time step. 
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The main idea of the UN02 scheme is to construct a non-oscillatory piecewise-parabolic 
interpolant Q{x) to a piecewise smooth function w{x) (see [46] for more details). On each 
segment containing the face Xi_^i e [xjjXi+i], the function Q{x) = is locally a 

quadratic polynomial and wherever w{x) is smooth we have 

Q(x) - w(x) = O(Ax^), ^(x±0) - ^ = O(Ax^). 

ax ax 

Also Q{x) should be non-oscillatory in the sense that the number of its local extrema does 
not exceed that of w{x). Since = Wi and = Wi+i, it can be written in 

the form 

X Xi ^ Xi^(^X Xi-^i) 



where d^^^Jyw) = Wi+i — Wi and D^_^iv is closely related to the second derivative of the 
interpolant since -Dj+i f = Ax'^q'f^i (x). The polynomial (x) is chosen to be one the least 

oscillatory between two candidates interpolating w{x) at (xj_i, Xj, x^+i) and (xj, Xj+i, Xj+2). 
This requirement leads to the following choice of D^^iv 

Di+iw := minmod(Aw, A+iw), DiW = Wi+i-2wi+Wi^i, A+iW = Wi+2-'2wi+i+Wi, 

and minmod(x, y) is the usual min mod function defined as: 

minmod(x,?/) = |[sign(x) + sign(?/) ] x min(|x|, ||/|). 

To achieve the second order (9(Ax^) accuracy, it is sufficient to consider piecewise linear 
reconstructions in each cell. Let L{x) denote this approximately reconstructed function 
which can be written in the form 

L(x) = Wi + Si (x — Xj) / Ax, Xj_i ^ X ^ Xj_^i. 

To make L{x) a non-oscillatory approximation, we use the parabolic interpolation Q{x) 
constructed below to estimate the slopes Si within each cell 



/ dQ 


dQ 




I dx 


' dx 


x=x+) 



Si = Ax X minmod | 

In other words, the solution is reconstructed on the cells while the solution gradient is 
estimated on the dual mesh as it is often performed in more modern schemes [7, 8]. A brief 
summary of the UN02 reconstruction can be also found in [32]. 

3.5. Time stepping. We rewrite the semi-discrete scheme (3.5) as a system of ODEs: 

dtw = C{w,t), w{0) = wq. 

In order to solve numerically the last system of equations, we apply the Bogacki-Shampine 
method [10]. It is a third-order Runge-Kutta scheme with four stages. It has an embedded 
second order method which is used to estimate the local error and, thus, to adapt the 
time step size. Moreover, the Bogacki-Shampine method enjoys the First Same As Last 
(FSAL) property so that it needs three function evaluations per step. This method is also 
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implemented in the ode23 function in Matlab [68]. A step of the Bogacki-Shampine 
method is given by 















h 


= 


) + |At„A;2,t„ + |At), 


^(n+l) 


= w(") + 


- At„ (|A;i + |A;2 + Ih) , 






f + Af 1 


W2 







Here w'-"^ ~ w{tn), At is the time step and w^2^^^ is a second order approximation to the 
solution w{tn+i)-, so the difference between and ■^2""''^^ gives an estimation of the 

local error. The FSAL property consists in the fact that ^4 is equal to ki in the next time 
step, thus saving one function evaluation. 

If the new time step At^+i is given by At^+i = p„A)f:„, then according to H211b digital 
filter approach [69, 70], the proportionality factor p„ is given by: 

where e„ is a local error estimation at time step t„ and constants /32 and a are defined 
as 

a = 1/4, /3i = /32 = l/4p. 
The parameter p is the order of the scheme (p = 3 in our case). 

Remark 8. The adaptive strategy (3.6) can he further improved if we smooth the factor 
Pn before computing the next time step Atn+i 

The function uj{p) is called the time step limiter and should be smooth, monotonically 
increasing and should satisfy the following conditions 

uj{0) < 1, w(+oo) > 1, uj{l) = l,uj'{l) = 1. 

One possible choice is suggested in [70]: 



u{p) = 1 + ftarctan 
In our computations the parameter k, is set to 1. 



(^) 
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Initial wavenumber n: 


1 


Gravity acceleration g: 


lms-2 


Final simulation time T: 


24 s 


Initial wave amplitude h\ 


0.2 m 


Undisturbed water depth d^: 


1 m 


Bathymetry oscillation amplitude a: 


0.1 m 


Low bathymetry oscillation wavelength ki. 


2m-i 


High bathymetry oscillation wavelength k2'- 


6 m^-^ 



Table 2. Values of various parameters used for the wave propagation over 
an oscillatory bottom test-case. 

4. Numerical results 

The numerical scheme presented above has already been validated in several studies even 
in the case of dispersive waves [32, 31]. Consequently, we do not present here the standard 
convergence tests which can be found in references cited above. In the present Section we 
show numerical results which illustrate some properties of modified Saint- Venant equations 
with respect to their classical counterpart. 

In the sequel we consider only one-dimensional case for simplicity. The physical domain 
will be also limited by wall-boundary conditions. Other types of boundary conditions 
obviously could also be considered. 

4.1. Wave propagation over oscillatory bottom. We begin the exposition of numer- 
ical results by presenting a simple test-case of a wave propagation over a static but highly 
oscillatory bottom. Let us consider a one-dimensional physical domain [—10, 10] which is 
discretized into N = 350 equal control volumes. The tolerance parameter S in the time 
stepping algorithm is chosen to be 10~^. The initial condition will be simply a bump 
localised near the center x = and posed on the free surface with initial zero velocity field 

Vo{x) = 6sech^(Kx), uo{x) = 0. 

The bottom is given analytically by the function 

d{x) = do + asin(A;x). 

In other words, the bathymetry function d{x) consists of uniform level c/q which is perturbed 
by uniform oscillations of amplitude a. Since the bathymetry is static, the governing 
equations (3.1) and (3.2) are simplified at some point. 

Hereafter, we fix two wavenumbers ki and k2 {ki < k2) and perform a comparison 
between numerical solutions to the classical and modified Saint- Venant equations. The 
main idea behind this comparison is to show the similarity between two solutions for mild 
bottoms and, correspondingly, to highlight the differences for stronger gradients. The 
values of various physical parameters used in numerical simulations are given in Table 2. 

Several snapshots of the free surface elevation during the wave propagation test-case are 
presented on Figures 3-9. The left image refers to the mild bottom gradient case (fci = 2) 
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Free surface elevation at t = 2.00 



Free surface elevation at t = 2.00 



(a) Low oscillations, ki (b) High oscillations, k2 

Figure 3. Wave propagation over an oscillatory bottom, t = 2s. 



mSV 
SV 




while the right image corresponds to the oscillating bottom (^2 = 6). Everywhere, the solid 
blue hne represents a solution to the modified Saint- Venant equations, while the dotted 
black line refers to the classical solution. 

Numerical results on left images indicate that both models give very similar results when 
bathymetry gradients are gentle. Two solutions are almost indistinguishable to grahic reso- 
lutions, especially at the beginning. However, some divergences are accumulated with time. 
At the end of the simulation some differences become to be visible to the graphic resolu- 
tion. On the other hand, numerical solutions on right images are substantially different 
from first instants of the wave propagation. In accordance with theoretical predictions (see 
Remark 6), the wave in modified Saint- Venant equations propagates with speed effectively 
reduced by bottom oscillations. This fact explains a certain lag between two numerical 
solutions in the highly oscillating case. We note that the wave shape is also different in 
classical and improved equations. 

Finally, on Figure 10 we show the evolution of the local time step during the simulation. 
It can be easily seen that the time adaptation algorithm very quickly finds the optimal value 
of the time step which is then maintained during the whole simulation. This observation 
is even more flagrant on the right image corresponding to the highly oscillating case. 



4.2. Wave generation by sudden bottom uplift. We continue to investigate various 
properties of the modified Saint- Venant equations. In this section, we present a simple 
test-case which involves the bottom motion. More precisely, we will investigate two cases 
of slow and fast uplifts of a portion of bottom. This simple situation has some important 
implications to tsunami genesis problems [43, 78, 29, 27, 50, 23, 57, 24]. 
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Free surface elevation at t = 5.00 



Free surface elevation at t = 5.00 




(a) Low oscillations, ki 



B 0.08 - 



"5^ 




(b) High oscillations, k2 



Figure 4. Wave propagation over an oscillatory bottom, t = 5s. 



Free surface elevation at t = 9.00 



Free surface elevation at t = 9.00 
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(a) Low oscillations, ki (b) High oscillations, 

Figure 5. Wave propagation over an oscillatory bottom, t = 9s. 

The physical domain and discretisation parameters are inherited from the last section. 
The bottom is given by the following function: 

1 2 



d{x,t) = do - aT(t)H(62 



X 



1 - e- 



-at 



where H(x) is the Heaviside step function [1], a is the deformation amplitude and h is 
the half-length of the uplifting sea floor area. The function T{t) provides us a complete 
information on the dynamics of the bottom motion. In tsunami wave literature, it is called 
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Free surface elevation at t = 16.00 



Free surface elevation at t = 16.00 
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(a) Low oscillations, k\ (b) High oscillations, ki 

Figure 6. Wave propagation over an oscillatory bottom, t = 16 s. 




Figure 7. Wave propagation over an oscillatory bottom, t = 18 s. 



a dynamic scenario [44, 43, 27, 50, 23]. Obviously, other choices of the time dependence 
are possible. Initially the free surface is undisturbed and the velocity field is taken to be 
identically zero. The values of various parameter are given in Table 3. 

Numerical results of the moving bottom test-case are shown on Figures 11-16. On all 
these images the blue solid line corresponds to the modified Saint- Venant equations, while 
the black dashed hne refers to its classical counterpart. The dash-dotted line shows the 
bottom profile which evolves in time as well. 
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Free surface elevation at t = 20.00 



Free surface elevation at t = 20.00 
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(a) Low oscillations, fci (b) High oscillations, k2 

Figure 8. Wave propagation over an oscillatory bottom, t = 20 s. 



Free surface elevation at t = 24.00 



Free surface elevation at t = 24.00 
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Figure 9. Wave propagation over an oscillatory bottom, t = 24 s. 

First, we present numerical results (see Figures 11-12) corresponding to a relatively slow 
uplift of a portion of the bottom {ai = 2.0). There is a very good agreement between two 
computations. We note that the amplitude of bottom deformation a/d = 0.25 is strong 
which explains some small discrepancies on Figure 12(a) between two models. This effect 
is rather due to the bottom shape than to its dynamic motion. 

Then we test the same situation but the bottom uplift is fast with the inverse character- 
istic time a2 = 12.0. In this case the differences between two models are very flagrant. As 
it can be seen on Figure 14, for example, the modified Saint- Venant equations give a wave 
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Time step size evolution 
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N , time step number 
(a) Low oscillations (ki) 



Time step size evolution 
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N , time step number 
(b) High oscillations (^2) 



Figure 10. Local time step size during the simulation of a wave propagation 
over an oscillatory bottom test case. 



Slow uplift parameter ai: 


2.0 


Fast uplift parameter a2- 


12.0 s"^ 


Gravity acceleration g: 


1 m 


Final simulation time T: 


5s 


Undisturbed water depth do: 


Im 


Deformation amplitude a: 


0.25 m 


Half-length of the uplift area b: 


2.5 m 



Table 3. Values of various parameters used for the wave generation by a 
moving bottom. 



with almost two times higher amplitude. Some differences in the wave shape persist even 
during the propagation (see Figure 16). This test-case clearly shows another advantage of 
the modified Saint- Venant equations in better representation of the vertical velocity field. 

On Figure 17 we show the evolution of the local time step adapted while solving the 
modified Saint- Venant equations with moving bottom (up to T = 5s). We can observe a 
behaviour very similar to the result presented above (see Figure 10) for the wave propaga- 
tion test-case. 

4.3. Application to tsunami waves. Tsunami waves continue to pose various difficult 
problems to scientists, engineers and local authorities. There is one question initially 
stemming from the Ph.D. thesis of C. Synolakis [74]. On page 85 of his manuscript, 
one can find a comparison between a theoretical (NSWE) and an experimental wavefront 
paths during a solitary wave runup onto a plane beach. In particular, his results show some 
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Free surface elevation at t = 0.50 



(a) t = 0.5 s 
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(b) t = l.Os 



Figure 11. Slow bottom uplift test-case (ai = 2). 
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(b) i = 5.0 s 



Figure 12. Slow bottom uplift test-case (ai = 2). 



discrepancy whose importance was not completely recognized until the wide availability 
of videos of the Tsunami Boxing Day 2004 [3, 59, 77]. In the same line of thinking, we 
quote here a recent review by Synolakis & Bernard (2006) [76] which contains a very 
interesting paragraph: 

"In a video taken near the Grand Mosque in Aceh, one can infer that the 
wavefront first moved at speeds less than 8 km h~^ , then accelerated to 35 
km h^^ . The same phenomenon is probably responsible for the mesmeriza- 
tion of victims during tsunami attacks, first noted in series of photographs 
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Free surface elevation at t = 0.50 



Free surface elevation at t = 0.90 
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(a) t = 0.5 s (b) i = 0.9 s 

Figure 13. Fast bottom uplift test-case (^2 = 12). 
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(a) t = 1.0 s 



(b) i = 1.5 s 



Figure 14. Fast bottom uplift test-case (02 = 12). 



of the 1946 Aleutian tsunami approaching Hilo, Hawaii, and noted again in 
countless photographs and videos from the 2004 megatsunami. The wave- 
front appears slow as it approaches the shoreline, leading to a sense of false 
security, it appears as if one can outrun it, but then the wavefront accelerates 
rapidly as the main disturbance arrives. " 

Since our model is able to take into account the local bottom roughness into the wave 
speed computation, we propose below a simple numerical setup which intends to shed some 
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Free surface elevation at t = 2.00 



Free surface elevation at t = 2.50 
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(a) t = 2.0 s (b) i = 2.5 s 

Figure 15. Fast bottom uplift test-case (^2 = 12). 
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Figure 16. Fast bottom uplift test-case (02 = 12). 



light on possible mechanism of the the reported hereinabove wave front propagation anom- 
alies. Consider a one- dimensional domain [—20, 20] with wall boundary conditions. This 
domain is discretised into = 4000 control volumes in order to resolve local bathymetry 
oscillations. The bottom has a uniform slope which is perturbed on the left side [x < 0) 
by fast oscillations which model the bottom roughness 



d{x) = cIq — xtan((5) + a[l — H(x) ] sin(A;x), 



(4.1) 
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Time step size evolution 
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Figure 17. Local time step size evolution during the simulation of a wave 
generation by moving bottom. 



Undisturbed water depth do: 


1 m 


Gravity acceleration g: 


1 m s^^ 


Bottom slope tan((5): 


0.02 


Oscillation amplitude a: 


0.1m 


Oscillation wavenumber k: 


20m-i 


Final simulation time T: 


19s 


Solitary wave amplitude A: 


0.3 m 


Solitary wave initial position xq'- 


-12.0 m 



Table 4. Values of various physical parameters used for the wave propaga- 
tion over a sloping bottom. 



where H(x) is the Heaviside function. The initial condition is a solitary wave moving 
rightwards as it was chosen in [74, 75]: 



Co rjoi^x) 



d{xo) + ?7o(x) 



3A _ v/6^(l + A) 7(1 + A) log(l + A)-A 

^ ~ ^'T+A' '° - V3T2A A • 

This configuration aims to model a wave transition from rough to gentle bottoms. The 
values of various physical parameters are given in Table 4. 

Then, the wave propagation and transformation over the sloping bottom (4.1) was com- 
puted using the classical and modified Saint-Venant equations. The wave front position 
was measured along this simulation and the computation result is presented on Figure 18. 
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The slope of these curves represents physically the wave front propagation speed. Recall 
also that the point x = corresponds to the transition between rough and gentle regions 
of the sloping beach. 

As one can expect, the classical model does not really see the rough region except 
from tiny oscillations. An observer situated on the beach, looking at the upcoming wave 
modelled by classical Saint- Venant equations, will not see any change in the wave celerity. 
More precisely, the slope of the black dashed curve on Figure 18 is rather constant up to 
the graphical resolution. On the other hand, one can see a drastic change in the wave front 
propagation speed predicted by the modified Saint- Venant equations when the bottom 
roughness disappears. 

The scenario we present in this section is only a first attempt to shed some light on 
the reported anomalies in tsunami waves arrival time on the beaches. For instance, a 
comprehensive study of P. Wessel (2009) [82] shows that the reported tsunami travel 
time is often exceeds slightly the values predicted by the classical shallow water theory 
(see, for example. Figures 5 and 6 in [82]). This fact supports indirectly our theory. 
Certainly this mechanism does not apply to laboratory experiments but it can be a good 
candidate to explain the wave front anomalies in natural environments. The mechanism 
we propose is only an element of explanation. Further investigations are needed to bring 
more validations to this approach. 

We underline that the computational results rely on sound physical modelling with- 
out any ad-hoc phenomenological terms in the governing equations. Only an accurate 
bathymetry description is required to take the full advantage of the modified Saint- Venant 
equations. 

5. Conclusions 

In this study we derived a novel non-hydrostatic non- dispersive model of shallow water 
type which takes into account large bathymetric variations. Previously some attempt was 
already made in the literature to derive shallow water systems for arbitrary slopes and 
curvature [22, 11, 49]. However, our study contains a certain number of new elements 
with respect to the existing state of the art. Namely, our derivation procedure relies 
on a generalized Lagrangian principle of the water wave problem [16]. Moreover, we do 
not introduce any small parameter and our approximation is made through the choice 
of a suitable constrained ansatz. Resulting governing equations have a simple form and 
physically sound structure. Another new element is the introduction of arbitrary bottom 
time variations. 

The proposed model is discretized with a finite volume scheme with adaptive time step- 
ping to capture the underlying complex dynamics. The performance of this scheme is then 
illustrated on several test cases. Some implications to tsunami wave modeling are also 
suggested at the end of this study. 

Among various perspectives we would like to underline the importance of a robust runup 
algorithm development using the current model. This research should shift forward the 
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accuracy and our comprehension of a water wave runup onto complex shores [62, 36, 37, 
30, 71, 33]. 
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